Statistics of Hartree-Fock Levels in Small Disordered Systems 
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We study the statistics of quasiparticle and quasihole levels in small interacting disordered 
systems within the Hartree-Fock approximation. The distribution of the inverse compressibility, 
given according to Koopmans' theorem by the distance between the two levels across the Fermi 
energy, evolves from a Wigner distribution in the non-interacting limit to a shifted Gaussian for 
strong interactions. On the other hand the nature of the distribution of spacings between neighboring 
levels on the same side of the Fermi energy (corresponding to energy differences between excited 
states of the system with one missing or one extra electron) is not affected by the interaction and 
follows Wigner-Dyson statistics. These results are derived analytically by isolating and solving 
the appropriate Hartree-Fock equations for the two levels. They are substantiated by numerical 
simulations of the full set of Hartree-Fock equations for a disordered quantum dot with Coulomb 
interactions. We find enhanced fluctuations of the inverse compressibility compared to the prediction 
of Random Matrix Theory, possibly due to localization of the wavefunctions around the edge of the 
dot. The distribution of the inverse compressibility calculated from the discrete second derivative 
with respect to the number of particles of the Hartree-Fock ground state energy deviates from 
the distribution of the level spacing across the Fermi energy. The two distributions have similar 
shapes but are shifted with respect to each other. The deviation increases with the strength of the 
interaction thus indicating the breakdown of Koopmans' theorem in the strongly interacting limit. 



I. INTRODUCTION. 



It is universally accepted that the principal signature 
of quantum chaos lifrthe statistics of Random Matrix The- 
ory (RMT), RefsB'El which is obeyed by the energy lev- 
els of chaotic systems. This is supported by semiclas- 
sical considerations] as well as many numericalEm and 
analytical^ examples. However realistic chaotic systems 
such as quantum dots, small metallic grains or the so 
called yrast levels in rotating nuclei involve interactions 
between many particles (electrons, nucleons, etc). Inter- 
esting problem then arises of how the interactions are 
expected to modify the RMT predictions. Several recent 
experimental and Jjhcpretical publications began to deal 
with this problemETO. 

Applications of RMT to noninteracting chaotic sys- 
tems, such as quantum dots, are concerned with sta- 
tistical properties of single particle quantities. Analo- 
gous and experimentally relevant in interacting systems 
arc quantities which characterize quasiparticles. In small 
disordered systems one can discuss the statistics of their 
energy levels, lifetimes and wave functions (respectively 
real and imaginary parts and the residues of the poles of 
the single particle Green's function). Experience gained 
in nuclear and atomic physics indicates that the Hartree- 
Fock (HF) method provides a very reasonable approxi- 
mate description of quasiparticle properties in finite sys- 
tems. The non degenerate HF particle-hole excitations 
form a convenient basis to describe low lying excitations 



of these systems. In this paper we adopt this descrip- 
tion and study statistical properties of the HF levels 
in small disordered systems. Properties of charged ex- 
citations have been probed experimentally by measur- 
ing Coulomb blockade in disordered quantum dots, cf. 
Refs.tnl The neutral particle-hole excitations can be 
measured tw studying acoustic phononlifl and microwave 
absorptionlL!]. 



II. THE HARTREE-FOCK APPROXIMATION IN 
WEAKLY DISORDERED SYSTEMS. 

Interacting electrons in a disordered system are de- 
scribed by the Hamiltonian 



H = h a palap + - ^ 

aj3 



+ t 



(i) 



where \a), \j3), . . . denote states of a single particle basis. 
The non interacting part of H is controlled by the one 
body Hamiltonian h a p representing the disordered sys- 
tem. In this work we are interested in the regime of disor- 
der for which the random matrix h a p can be viewed as de- 
scribed by the rules of RMT. The interaction V a /3-yS is not 
random and in a quantum dot, for example, represents 
matrix elements of the Coulomb or screened Coulomb in- 
teraction. In our discussion however we will regard it as 
a general matrix and will be able to draw conclusions for 
wide classes of possible Va^s- 
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Our main approximation will be to treat the Hamil- 
tonian (|l|) in the HF approximation. The HF equations 
are 

X>g^*(/3) = £ilM«) , (2) 
P 

with 

h% F) =h a0 + Yl V <L>wPP><*' . (3) 

a'0' 

or symbolically h^ HF ^ = h + tr(V A p). Here V A is the an- 
tisymmetrized interaction matrix V A p S = V a pys — V a fisj 
and the trace is taken over the second pair of the indices 
of V A . The self-consistent density matrix p is given by 



Pa0 = y, Maw h w)- (4) 

h^holes 

Henceforth we will denote by h, h' andp, p', etc., the hole 
(occupied) and particle (unoccupied) levels respectively. 

It is common to view the HF approximation as aris- 
ing from the variational minimization of the total energy 
of the system. In this interpretation the direct physical 
meaning of the energies e, and the wave functions ipi ( a ) 
remains obscure and-pne must use the so called Koop- 
mans theorem, Ref£3. We recall however that the HF 
equations can also be derived from an approximation to 
the equation of motion for the one particle Green's func- 
tion 



G(a, /3; u>) = —i 



dt 



($ (N)\Ta a (t) a+(0)|$ o (iV)) 



i 

£ 



($ (AQ| aa |$. t (jV + 1))(MN + l)\ap$o(N)) 
Ei(N + 1) - E (N) - lu - iO 

(^o(N)\a+\^(N - l))($i{N - l)|a a |$ (JV)) 
E (N) — Ei(N — l)-u + iO 



(5) 



The HF energies e p and Sh and the corresponding wave functions ip p (a) and ij)h( a ) are then approximate energies and 
wave functions of respectively quasiparticles and quasiholes, 



e p ~ E P (N + 1) - E (N) 
xP p {a) ~ {^{N)\a a \%(N + 1)) 



s h ~ E (N) - E h (N - 1) 
ip h (a) ~ ($ h (N - l)\a a \$ (N)) . 



We wish to study the statistical properties of the set 
£i and the wave functions ipi(a) which follow from the 
random nature of h a p. Ideally one would like to be able, 
starting from the probability distribution P{h), to deter- 
mine the joint probability distribution P(e\, £2, £n, ■■■) 
and similar distribution for ipi (/?) and on its basis to pre- 
dict various correlation properties of these quantities. In 
the present paper we address a much simpler problem of 
the repulsion pattern of neighboring pairs of Si and its 
application to the addition spectra of quantum dots. We 
will present simple analytic approximations and perform 
numerical investigations to check their validity. 



(6) 



This interaction allows for a trivial exact solution which 
is reproduced by the HF equations. They are solved by 
the eigenfunctions ipi ( a ) °f the one body Hamiltonian 

h a ^f\(3)=ef ) i>f\a), (7) 
and have the following eigenvalues 



.(0) 



+ V (N-1) 



VnN 



Til = 1 

m = 



(8) 



III. THEORETICAL CONSIDERATIONS. 

A. The Constant Part of the Interaction. 

The simplest limiting case is a constant interaction 
V(\r — r'|) = Vq which will serve as a reference point for 
our discussion. In the following we consider the spinlcss 
case. One has for the matrix elements in (nu 



where N is the number of particles and rii is the occupa- 
tion number of the i'th level. It is clear that the lowest 
energy solution is obtained by occupying the lowest N 
non-interacting states. The gap which separates the en- 
ergies of the occupied (hole) and empty (particle) levels 
for Vq > results from the absence of the contribution 
from the exchange term in the latter. In a sense one can 
say that the empty levels do not create the exchange hole 
which decreases (increases) the level energy for positive 
(negative) Vq. 
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From the known distribution of 's it is easy to cal- 
culate the statistical properties of the HF levels e , . Since 
the spacings of neighboring levels are s = Ae^ for lev- 
els below and above ej their distribution is given by 
the ordinary Wigner surmise (e.g. P(s) = P\y(s) ^ 
s exp[— 7r( 2< s s> ) 2 ] f° r GOE). The spacing between the 
levels lying across the Fermi energy is s = Vo + Ae^ ' 
and therefore its distribution is given by the Wigner sur- 
mise shifted by the value Vo, 



P(s)=P w (s-V ) 



(9) 



The gap in the distribution of the level spacings across 
the Fermi energy is reflected in the differences of the to- 
tal ground state energy Eq of the system as a function of 
the particle number 



E (N) 



N 



.(0) 



Vn 

Y N ( N ' 

-(0) 



1) 



E (N + 1) - E Q (N) = e%> +1 + V N = e N+1 , 
D 2 {N) = E (N + 1) - 2E {N) + E {N - 1) 



-(0) 
~'N+ 



Vo — £n+i — £n 



(10) 
(11) 
(12) 



Here we defined the quantity D 2 (N) which is essentially 
the inverse compressibility of the dot. Its distribution 
has been measured in transport experiment, of quantum 
dots in the Coulomb blockade regime, Refs.BTj. Eq. ( O ) 
is Koopmans theorem which is exact in this case. It 
leads directly to the expression ( |l2| ) of D<x in terms of 
the energy difference between the two HF levels across 
the Fermi energy, a fact which motivates the study of 
this spacing. The Coulomb gap in D 2 does not fluctu- 
ate in the limit of constant interaction. For more realis- 
tic interactions the Coulomb blockade gap must undergo 
fluctuations in addition to the fluctuations of the single 
particle energies Si. Study of these fluctuations will be 
one of our principal goals. 



B. Realistic Interaction — Repulsion of the 
Hartree-Fock Levels. 



For small s (s <C A - the mean level spacing) the 
well known behavior of the distribution P(s) for random 
Hamiltonians (i.e. linear for GOE, quadratic for GUE, 
etc.) can be derived by considering the repulsion of close 
pairs of levels. In this subsection we generalize this anal- 
ysis to the case of two close HF levels. The new element 
in the HF problem, Eq. (|J), is the presence of the non 
linear self-consistent term X^a'fl' ^ pp'PP'a' which im- 
plicitly depends on the realization of the random part 

Let us assume that for some realization of h a p the HF 
Hamiltonian h!fL F ^ has two closely lying levels, £2 > £1, 
so that Ae = £ 2 — £1 <C A. We wish to investigate 



how Ae reacts when the random h is varied by Sh with 
A > \Shu\ ~ Ae for I, J= 1,2. 

The variation of h causes a change Sh^ HF ^ = Sh + 
tr(V A 8p) in the HF Hamiltonian where the second term 
is due to the induced change of the self-consistent density. 
Without this term one would get for the new spacing the 
standard result s = 2^(Ae + 5h 22 ~ Shu) 2 + 4|<5/ii 2 | 2 , 
which for future reference we shall rewrite in the form 



B\ 



R2 

-°3 



2IBI 



(13) 



The vector B is a projection on the Pauli matrices of the 
2x2 Hamiltonian (h + Sh)u in the subspace of the two 
close levels, I, J — 1, 2, 

B = — tr[<r(h + Sh)] , cr, — the Pauli matrices. (14) 

The level repulsion at small spacings is a consequence of 
the proportionality of s to the square root of the sum of 
squares of real quantities and can be deduced by evalu- 
ating, for small s 



P(s) = J 8(s~ «(B))P(B)dB . 



(15) 



The familiar behavior P(s) ~ s for GOE and P(s) ~ s 2 
for GUE (s < A) follows then directly for any P(B) 
which does not vanish for small B = |B|. The (unjusti- 
fied) evaluation of (|I1f ) for all values of s taking P(B) to 
be of a Gaussian form results in the Wigner distribution 
for the case of real h. 

In the Appendix we discuss how the result ( |l3| ) is mod- 
ified by the presence of the self-consistent t erm tr (V Sp). 
Under the conditions discussed in Section III C we show 
there that : 

1. When the two close levels are both occupied or both 
empty the distribution of their level spacings for small 
s <C A is ~ s or ~ s 2 as in the non interacting case. 

2. When the two levels are on opposite sides of the 
Fermi level, i.e. one occupied and one empty the expres- 
sion (Till) changes to 



B J m(B, J 



(16) 



Here we use the dyadic notation [J -m] a = J2b=i Jabtnb- 
As in the non-interacting case the "vector" B contains 
the information about the random part h + Sh of the 

HF Hamiltonian whereas the matrix J is the projec- 
tion on the Pauli matrices of the interaction matrix 
Vijkl, (I,J,K,L = 1,2) in the subspace of the two 
close levels tpi and ip% (cf. Appendix for the precise def- 
inition of V, B and J). The vector m(B, J) is a unit 
vector |m| = 1 and is a solution of the self-consistent HF 
equations for the two levels 



m x (B— J ■ m) = . 

In this notation the total HF energy is 



(17) 
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£ = B m — -m J m . 



(18) 



The full investigation of the solutions of the above 
equations and the expression is found in the Ap- 
pendix. Here we will concentrate on the simplest case 
when the matrix of interactions J a b is degenerate, i.e. 
has equal eigenvalues. As is shown in the Appendix such 
a degenerate J a f, matrix corresponds to a complete ab- 
sence of degeneracy (such as spin) of the close HF levels. 
We will furthermore restrict ourselves to real Hamilto- 
nians h + 5h. Then the vector B has only two com- 
ponents B± and B3 and one needs to consider only the 
corresponding components of J a b with a,b = 1,3. For a 

degenerate J a b the vector J m is parallel to m and the 
two level HF equation ([l7]) becomes a simple linear rela- 
tion miBj, — m^Bi = without any dependence on the 
interaction. The level spacing (16) however still contains 
the interaction term. Using the normalization condition 



771.3 = 1 one obtains 



s = 2\B + J\ 



(19) 



where B = \J B\ + i?| and J is the degenerate eigen- 
value. The explicit expression for J is 



J = {Vr. 



V122O/2 



(20) 



For a positive J which corresponds to the case of a 
repulsive interaction the equation for the spacing ( |l9| ) 
describes a conical surface separated by a distance 2J 
from the {B\,B^) plane. This means that the probabil- 
ity for a given s has a gap of the magnitude 2J above 
which it must rise linearly 



P(s;J) 





(s - 2J)f(s) 



s < 2 J 
s > 2 J 



(21) 



where /(s) is an undefined function which is finite at 
s = 2J. For general reasons one should expect that f(s) 
vanishes at s — > 00. 



The result (21) implies that the distribution P(s; J) is 
qualitatively the same as the shifted Wigner distribution 



2Jt(s - 2J)e~7( — ) S >2J v ' 

similar to what was found in the case of the constant in- 
teraction, Eq. (||). There is however a crucial difference 
which we will now try to elucidate. 



As is seen from Eq. ( p0| ) the degenerate eigenvalue J is 
invariant under unitary transformations in the subspace 
spanned by the two given close states ipi and ip2 ■ As long 
as one examines realizations of the random h for which 
these two particular states stay close in energy they mix 
strongly only between themselves and J stays the same. 
The distribution of the level spacings for such selection 
of h's is given by ( |2l| ) with the fixed J as in the constant 
interaction case. However variations of h may bring an- 
other pair of levels, say, ipi and -03 at close distance on 
both sides of et. For such levels the value of J may 
be completely different and given by the corresponding 
matrix elements (V1313 — Vi33i)/2. In calculations of the 
overall P(s) for small s one should therefore average over 
the distribution P{ J) of all J's 



P(s) = J P(J)P(s;J)dJ 



(23) 



For a constant interaction J = Vo for any pair of 
states tpi, ipj and therefore P(J) is a delta function cen- 
tered at Vo. For an extreme short range interaction 
V(\r — r'|) = VoS(r — r') all values of J are zero and 
do not fluctuate. One expects therefore that for a gen- 
eral interaction the average value of J is coming mainly 
from the very long range component while its fluctuations 
reflect the middle range components of V. 

The main feature of matrix elements (jfy of a constant 
interaction is that they have the same form in any (or- 
thonormal) basis of the single particle states ?p a which are 
used to calculate Va^s- In other words for a constant in- 
teraction the matrix elements Vap^s are invariant under 
the unitary group U(M) of all transformations in the en- 
tire single particle Hilbert space of the problem. Here M 
is the dimensionality of this space (typically M — > 00). 
This invariance of Vapyg for a constant interaction is the 
reason that P( J) is a delta function in this case. In or- 
der to develop a theory of P(J) for a general interaction 
one must understand the statistical properties of the sin- 
gle particle wave functions-solutions of the HF problem. 
The matrix elements in (|20| ) are defined with respect to 
such wave functions. Although we do not have an an- 
alytic theory of P{J) we will present below numerical 
investigation of this function and the validity of the ex- 
pression of the type (p3|). Numerically we find that the 
distribution of J is approximately Gaussian. Adopting 
this finding into expression ( |23] ) and assuming the form 
(^||) for P(s; J) one obtains 



P(s) = (5 



a 



a + P 



0?(s-2Jo) 2 



■ erfc 



a + /3 Va+~P 



(s - 2 J ) 



T (s-2J ) 2 



(24) 
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where a — l/8aj 2 and (3 — 7r/4A with A being the 
mean level spacing at the vicinity of the Fermi energy. 
Jo and a j are the mean and standard deviation of P(J) 
respectively. 



C. Validity of the Two Level Treatment - Why Only 
the Statistics of Spacings Across the Fermi Level is 
Effected by the Interaction. 

The above results are valid as long as the spacing be- 
tween the given two levels is much smaller than the dis- 
tance to other levels. This restriction is needed in order 
to be able to isolate the two levels from the rest. For re- 
pulsive interactions the distance between the pair of lev- 
els on different sides of the Fermi energy is determined by 
the size of the matrix element J (cf.,Eq. [l9]). Therefore 
the condition of the validity of our treatment for such 
levels is J <C A. As we will see below numerical evidence 
indicates that at least qualitatively the expression ( p3| ) 
remains correct also for a much larger J ~ A. For an 
even stronger interaction the two level treatment ceases 
to be valid and one must account for the reaction of dis- 
tant levels to the changes of the self-consistent potential 
tr(V A Sp). 

Here we wish to add the following remark. For a 
constant interaction our result is trivially valid for any 
strength. On the other hand for any given interac- 
tion one can extract a constant part, i.e. V a p-ys — 
VoSa-fSps + Uap^s, where Ua^s = V a p~fS — V 5 aj 5ps- 
Since the HF wavefunctions are independent of Vq one 
may try to solve the problem using U first and then add 
Vq . This procedure is not unambiguous and must be dic- 
tated by the physics of the problem. If P(J) vanishes for 
J < J c , i-c in the presence of a hard gap (like in a quan- 
tum dot with Coulomb interaction), it is natural to take 



Vq = J c . The condition of validity is then J r <C A where 
J r is the typical gap due to the residual interaction U. 

When there is no minimal value for J any subtraction 
is bound to produce U which has both attractive and re- 
pulsive components. This fact may introduce fundamen- 
tal differences between the solutions to the HF problem 
obtained using the interaction V and the one derived in 
the manner indicated above. Most notably using the lat- 
ter procedure one will find cases for which there exists 
enhanced probability for the two levels across e/ to be 
close to each other [or at distance Vq after the addition 
of the constant part, see result ( Al5| ) of the appendix]. 
However if the width oj of P(J) is much smaller then its 
mean J the use of Vb = Jo — &j will cause U to change 
sign only in a small number of cases. Consequently we 
may use in such a case uj < A as our criterion for the 
validity of the two level treatment. 

For pairs of neighboring levels which lie on the same 
side of the Fermi energy the change in the density ma- 
trix Sp vanishes as long as only the two level subspace 
is considered. Consequently there is no self-consistent 
term V A 6p present and one recovers the results of the 
non-interacting case. There are two types of corrections 
to this two level treatment which must be considered. 
One is the non-zero contribution of the distant levels to 
dp. Another is the correction to the energies of the two 
close levels due to virtual transition to distant levels. 

To estimate the first correction we use the first order 
result 



Sp a f3 = ^2 

ph 



(p\6h + V A Sp\h) 



£h 



(25) 



which can be solved to obtain the RPA expression 



6p a = ^2 



OL f3 



Thus we find 



6p. 




dh 



(p'\V A \h 



h.c. 



V 



(27) 



The contribution to 
but 6h( HF ^ ~ 5h[l + 



where here V ~ J2 P ',h' V PP'hh' 
$h( HF ^ is therefore not just 5h 
V/ (V + A)]. However this still means that for such levels 
5h( HF ^ ~ Sh for any strength of interaction V, the sole 
role of which is to renormalize the random part 8h. 

The corrections to the energies of the close levels due 
to transitions to distant levels are 



E 

ph 



( P \Sh\h) 



(26) 



E T =£1 



E 

i/1,2 



(HF), 2 



£l - £j 



1,2 



(28) 



where the e/'s are the energies obtained in the two level 
treatment. But we have just shown that 5h^ HF ^ ~ 5h. 
Therefore the correction term in expression ( p8| ) is of or- 
der (5h) 2 /A <C Sh and can be disregarded for any V. 
Thus we expect that the two level treatment for a couple 
of occupied or empty levels gives correctly the small s 
behavior of P(s) for any interaction strength. To stress, 
the difference between this case and the case of two levels 
on both sides of e/ is that there 5p also included a " non 
perturbative" zero order term coming from the mixing of 
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the wave functions of the occupied and the empty states. 



D. Addition Spectrum vs Excitation Spectrum. 

The HF energies are interpreted as energies of quasi- 
particles (for e% > £/) or quasiholes (e* < £/). They are 
excitations of a system with one added or one subtracted 
particle, i.e. E^N ± 1) - E {N) = ±£j. The excita- 
tions with the same particle number Ei(N) — Eq(N) are 
described in the HF approximation by solutions repre- 
senting determinants with the same particle number N 
which are orthogonal to the HF ground state determi- 
nant. In the case of a constant interaction such exci- 
tations are simply particle-hole excitations of the non- 
interacting problem. The distribution of their spacings 
does not show any Fermi energy related gap and coincides 
with the non shifted Wigner distribution. This is also 
true in the two level HF model, Eq. ( |l7| ) . Let us demon- 
strate this in the simple case of the degenerate J . The 
HF equation ( |l7| ) in this case is mi -B3 — m^Bi — which 
together with the normalization condition m\ + m| = 1 
produce two solutions with different total HF energies, 
Eq. (plj|), £ = ±B — J/2. These solutions represent 
the ground and excited states of this model which have 
the same number of particle. The linear dependence of 
the difference A£ = 2B indicates that at least as long 
as the two level treatment of the level repulsion is valid 
the HF energy spacing distribution between such states 
obeys P(A£) ~ A£ for small A£ without any gap, as in 
the non-interacting systems. 

Yet another way to obtain this result is to consider 
the spacings between neighboring particle or hole levels. 
They correspond, within the HF approximation, to the 
distances between adjacent excited states of the system 
with N + 1 and N — 1 electrons respectively. By the ar- 
guments given in the preceding subsection these spacings 
follow the Wigner-Dyson statistics. 



E. A Schematic Model — Keeping Only the Average 
Interaction Matrix Elements. 



For a given realization of the random h a p let us con- 
sider the Hamiltonian in the eigenbasis of h, for which 
h a p — Sa^d a /}. In this random basis also the matrix ele- 
ments of the interaction V a p-ys are r^mdom. Their statis- 
tical properties are known, cf., Refs.E3 and are as follows. 
Only the matrix elements V a p a p and V a pp a have non 
zero averages. Their distributions are narrow with the 
width behaving like 1/M in the Random Matrix Theory 
(M - the size of the single particle space) and like 1/g 
in the random potential theory (g - dimensionless con- 
ductance). Based on these properties one is tempted to 
approximate the interaction by retaining only the matrix 
elements with non zero averages, i.e. to assume that 



Vcfrs = ScrySps V} ap} + 5 a5 5p 7 V^ al3) , (29) 

with V^ ^ = V2 . Such a model has an easy exact 
solution. The Hamiltonian 



H = ^2 £ a } n a + - £ J a p n a hp , 



(30) 



a,P 



where J a p = v} 01 ^ — V 2 ^ Q/3 \ has exact eigenstates given 
by the eigenfunctions of the occupation operators h a 



ijj = \ni,n 2 , ■ ■ ■ ,n k , ■ ■ •) , 
with the corresponding eigenenergies 



(31) 



(32) 



Like in the case of the constant interaction the exact re- 
sults are reproduced by the HF approximation. The HF 
equations are 



4 0) +E J ^^ 



i a ) - £ Ja(3 PaP <j>i{fi) = £i4*i{u) 

(33) 



and are solved by <pi (a) = 5i a . Therefore 

J0) -£•/,,, (34) 



£i = £, : 



where the sum is over the occupied levels. 

We are interested in the ground state ofpthe Hamilto- 
nian (]3^). A general argument due to Ref.Eil guarantees 
that at least for positive definite (repulsive) interactions, 
which we assume below, the HF ground state must be 
comprised of the N lowest energy single particle HF lev- 
els (j>i, ■ ■ ■ , (f)^. While the ground state of the constant 
interaction model is obtained by filling the N lowest non- 
interacting states this need not be the case for the present 
model. Consider 



_ (o) (0) , T 

£JV+1 — £N — £ N+1 — £ N ' JN+l.N 
JV-1 

+ £ (JN+l,k — JN,k) 
k=l 



(35) 



Although £^ +1 — + Jn+i,n is positive definite the 
sum in ([35]) is over 2A^ — 2 random variables and can be 
large and negative. Consequently it may happen that 
£n+i < £n in variance with the above condition on the 
ground state. In such different occupation pat- 

tern must be sought (we note that even if the condition 
is not violated there may exist other solutions that are 
consistent with it and which have lower energies). We 
expect such crossings of levels across the Fermi energy to 



.(0) 



G 



take place when A + Jo — V^Ncrj, where Jo and aj are 
the mean and standard deviation of J a p respectively. 

Similar crossings may occur for levels below or above 
e f when A ~ y/2N aj since for them 



iV 



-l - 



— e 



(0) 
3-1 



Y^( J Jk - Jj-i,k) ■ (36) 



fc=i 



This shuffling of levels tends to reduce the correlation 
between the energies of neighboring states. In particular 
one expects to find weaker level repulsion when either N 
or the interaction strength (<tj) are increased. 
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FIG. 1. Probability density distributions of the spacings 
between a) the last two hole levels below £/, b) the two lev- 
els across ej and c) the first two particle levels above e/. 
The results were derived by solving the HF equations of the 
schematic model Eq.(|30|). and J a/ 3 were generated by the 
random potential model with Coulomb interaction, described 
in Section for a dot with 15 electrons, W/t — 1.2 and 
U/t — 1.2. We also included a constant interaction part of 
strength Vo = 6Ao. The spacings are measured in terms of 
the non-interacting mean level spacing Ao. 

To verify the above discussion we have calculated the 
HF energies (34) using numerical values for and J a p 
generated by the random potential model described in 
the next Section. We took care to choose the lowest en- 
ergy solution. The results for P(s) below, across and 



above Sf are shown in Fig. [I]. These distributions dif- 
fer significantly from the exact HF distributions calcu- 
lated while retaining the off-diagonal elements of V, cf., 
Figs. The most prominent incorrect feature is 

the absence of level repulsion in P(s) for the levels be- 
low or above £/. One can attempt to correct this fea- 
ture by using the average rather than the exact matrix 
elements which enter J a p. .As it follows from the ran- 
dom potential model, Refs.E3, such average matrix ele- 
ments are functions of the corresponding eigenenergies, 
i.e. J Q/ 3 = j(\s a — ep\) with j(x) a known function which 
is approximately — In x in two dimensions. Although this 
model reproduces the density of states in the vicinity of 
the Fermi energy it is not expected to account for the 
correct correlations between neighboring levels as mani- 
fested in P(s). 



IV. COMPARISON WITH NUMERICAL 
RESULTS. 

In order to substantiate the results of our analytical 
two level treatment we have numerically solved the com- 
plete set of HF equations (Q) derived from a tight-binding 
Hamiltonian for a disordered two dimensional quantum 
dot. With the labeling of the sites by a double index 
the one-body Hamiltonian is given by 



i,3 



i,j u h3 



i,3 



R +1 



j">i,3 



h3 



+ h.c. 



(37) 



where £jj is the energy of the site and t is a con- 
stant hopping matrix element. Each of the energies £jj 
is chosen randomly from a Gaussian distribution with 
the standard deviation W/2. We assume repulsive 1/r 
interaction 



V = U 



\ - 4,j a k,i a k,iai,j 



i,j,k,l 



r k,l\/ b 



(38) 



where b is the lattice constant and U — e 2 /b. 

The dot was approximated by a grid of M = 20 x 20 = 
400 sites with hard wall boundary conditions. Most of 
the numerical data was obtained for a dot filled with 
N = 15 — 17 spinless electrons and a disorder strength 
W = 1.2 1. Under such conditions the dot was in the 
diffusive regime and the levels in the vicinity of £/ ex- 
hibited RMT statistics in the non- interacting limit. For 
the low filling that we used the energy band was approx- 
imately parabolic and the Fermi energy for a clean non- 
interacting system was £/ = AntN/M. A convenient di- 
mensionless measure of the strength of the interaction is 
r s = (e 2 /a)/ef where a = y/ M/irNb is the average inter- 
particle distance. In our case r s = (U/t) y/M/lQitN ~ 
0.7 U/t. Below we describe our results for U/t = 0.2—1.6. 
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Henceforth energy is quoted in units of the observed non- 
interacting mean level spacing at the Fermi energy Ao 
which was found to be larger by 7% then the clean value 
Airt/M. 




FIG. 2. Probability distributions F(s) = / Q a P(s') ds' of 
spacings between the last two hole levels just below Sf. The 
heavy lines depict from top to bottom the numerical results 
for U/t = 0, 0.4, 0.8, 1.2 and 1.6. The solid curves present the 
best fit to the data assuming a Wigner function with renor- 
malized mean level spacing A. The inset contains P(s) for 
U/t = 1.6 (histogram) together with its best fit to a Wigner 
function. 




FIG. 3. Same as Fig|2j but for the spacing between the first 
two particle levels just above the Fermi energy. 

The distribution functions F(s) = P(s') ds' of the 
spacings between the last two occupied (hole) levels and 
between the first two empty (particle) levels are shown 
in Figs. H and [| They are compared with the best fit 
to an integrated Wigner function with a renormalized 
mean level spacing A. These and the following results 
were obtained by averaging the HF spectrum of 15,16 
and 17 electrons over 450-500 realizations of the disorder. 



While the distributions vanish quadratically for small 
spacings increasing deviations from the Wigner function 
are observed when the interaction becomes stronger. We 
find enhanced probability for the occurrence of spacings 
smaller and much larger then A for large values of U/t. 
The renormalized mean level spacing A also increases 
with the strength of the interaction. The width of the 
fitted Wigner function <j(Pw) — 0.52 A is presented in 
Fig. It grows approximately linearly with U ft. The 
mean level spacing between adjacent levels further away 
from the Fermi energy decreases with the distance from 
e/ and approaches 1. 

The distributions of the quantity J defined in Eq. ( |20| ) 
are depicted in Fig. ^[ They are approximately Gaussian 
for small values of the interaction strength but develop 
asymmetry towards the high J end when U/t is increased. 
The mean of the distribution scales with the interaction 
as (J) ~ 1.7 U/t. For the width a (J) = ^(J 2 ) - (J) 2 
we find ct(J) ~ 0.16 U/t + 0.13 {U/t) 2 over the range of 
parameters studied (see Fig. |7|). These fluctuations are 
responsible for the smearing of the distribution of spac- 
ings across e/ as will be shown below. A lengthy calcula- 
tion using the random vector model (RVM) gives for our 
system (J) = 2.1 U/t and a (J) = 0.032 U/t. The fact 
that the RVM result for (J) is larger then the observed 
one reflects the tendency of the system to prefer a non- 
uniform density distribution that reduces the Coulomb 
energy. We believe that this is also the reason for at least 
part of the enhancement of the actual fluctuations rela- 
tive to the RVM predictions. Typically we found the HF 
eigenfunctions (both particles and holes) to have large 
amplitude along the periphery of the dot, as expected 
from simple electrostatic considerations. The localizing 
effect of the interaction is evident from Fig. || where we 
present the distribution of the inverse participation ratio 
/ = Alb 4 J2ij ip 4 {h3) averaged over the 10 particle and 
10 hole levels around e/. 
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FIG. 4. Probability density distributions of the parameter 
J for various interaction strengths. 
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FIG. 5. Probability density distributions of the inverse par- 
ticipation ratio of the 10 hole and 10 particle levels around Sf 
for various interaction strengths. 
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ing to its definition (|lj) and using the HF many body 
ground state energies for 15,16 and 17 electrons. We also 
compare them to the analytic estimate ( p4| ) . For its eval- 
uation we used A that interpolates between the values 
found from the Wigner functions fitting P(s) below and 
above £/ and the numerical results for (J) and tr(J). The 
two distributions evolve from shifted Wigner functions 
at small values of U jt to an approximate Gaussian dis- 
tributions as the interaction strength is increased. The 
crossover occurs when r s ~ XJ It ~ 1 . Around this point 
the fluctuations of the Coulomb gap a(2J) are compara- 
ble to the width of the Wigner function describing P(s) 
at the vicinity of the Fermi energy (see Fig. Conse- 
quently the latter is smeared into a new more symmetric 
and broader distribution. It is evident from Fig. |^ that 
Koopmans' theorem breaks down for strong interactions 
i.e. P(Ae) ^ P(Da). However it seems that the two dis- 
tributions may be viewed as shifted versions of each other 
having similar shapes but somewhat different widths (see 
Fig. ^). The shift of P(D 2 ) towards lower values is ex- 
pected and is due to the change of the occupied levels 
in response to the additional electron. This rearrange- 
ments, which is neglected by Koopmans' theorem, tends 
to lower the electrostatic charging energy. For reasons 
that are not clear to us our analytic estimate for P(Ae) 
fits rather well P(D 2 ). 
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FIG. 7. The standard deviations of 2J, Pw - the Wigner 
distribution that hts best the distribution of spacings in the 
vicinity of e/, D2(16) and Ae - the spacing between the two 
levels across the Fermi energy. The inset depicts the normal- 
ized fluctuations oiD-z)! {D2) (full diamonds) and <r(Ae)/(Ae) 
(empty diamonds). The results are for a dot with 15-17 elec- 
trons and disorder strength W = 1.2 1. 



FIG. 6. Probability density distributions of the spacing be- 
tween the levels across e/ (histograms) and of D2 (16) (broken 
lines) . The solid curves correspond to the estimate (M . 

In Fig. H we present the distribution for the spacing 
Ae between the HF levels across the Fermi energy to- 
gether with the distribution of Z? 2 (16) calculated accord- 



We repeated the numerical calculations for the same 
dot but with stronger disorder W/2 = 1.6 1 and fewer 
(10-12) electrons. All of the effects reported above have 
been observed for this case as well. The deviation of P(s) 
above or below Sf from the Wigner function were more 
pronounced. The fluctuations of J also increased and 
we found o-(J) ~ 0.25 U/t + 0.06(U/t) 2 with enhanced 
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asymmetry in the shape of the distribution for strong in- 
teractions. The discrepancy between P(Ae) and P{D 2 ) 
at large r s persisted although the shapes of the two distri- 
butions and particularly their width were closer for this 
dot then for the one described above. Fig. || summa- 
rizes the dependence of the fluctuations of the different 
quantities on the strength of the interaction. 
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FIG. 8. Same as Fig. | but for a dot filled with 10-12 
electrons and disorder strength W = 1.6 1. 



Our calculations were done for a fixed number of parti- 
cles. In order to facilitate comparison with a fixed chemi- 
cal potential ensemble the Fermi energy of each spectrum 
in the ensemble was shifted to zero. The resulting density 
of states is plotted in Fig 0. 



U/t=0.2 




FIG. 9. The density of states near e/. The results are 
for a dot filled with 15-17 electrons and disorder strength 
W = 1.21 



APPENDIX A: THE TWO LEVEL HARTREE-FOCK PROBLEM 



In subsection [II C we outlined the way in which the density matrix p changes under a variation Sh of the random 
part of the HF Hamiltonian. We argued that while in the case of two close levels that are both occupied or empty 
the term tr(V A Sp) does not alter the small s behavior of P(s) it must be included in a self-consistent manner for a 
couple of close level residing on both sides of £f. In this appendix we will concentrate on the latter scenario assuming 

Let us denote by ipi the eigenfunctions of h^ HF ^ and by fa the eigenfunctions of h^ HF ^ + 5h^ HF \ Then in the basis 
of 4>i and ip2 the HF problem for two close levels on both sides of Ef is 



4 0) + Shu + J2 VmjSpji ) a + 6h 12 + J2 Vii2j5p.ji )b = ea , 

J,J=1,2 J \ 1,7=1,2 J 

Sh *2 + E VznjSpJi ) a + 4 0) + Sh 22 + E V 2 i2jSpji )b = eb , 

/,,/=!, 2 / \ /,J=1,2 / 



(Al) 



with Spji the 2x2 matrix 

,._^^+ ./../.+ _/"!« 

a*b \b\ 



sp = Mt-Mt = [ |a| " 1 S ) . ( A2 ) 
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and the normalization condition |a| 2 + |6| 2 = 1. The relation between the matrix Vijkl, I — 1,2, etc., and the original 
Vap-yS will b e di scussed below. We will also need the part of the HF energy the minimization of which, with respect 
to 0i , gives dAlj) 



E HF [4n] = tr \ h + tr(V P ) + Sh 



n<Pi 



-tr 



(A3) 



Transformation to an effective magnetic problem. The solution of the set of algebraic equations (Al) is 
particularly simple if one uses the expansion in terms of Pauli matrices, 



1 



Mi = o(l+ m " CT ) 



1 



(1 + m • tr) 



(A4) 



with the "pure state" conditions |m| = |mo| = 1. Inserting this into (A3) and transferring to the left all terms which 
do not depend on m one obtains 



£ = B ■ m -m- J m 



(A5) 



The notation in ( |A5| ) stands for 



£ = E HF - -tr 



B 



Jab 



h + tr(Vp) + Sh +-tr{tr V(m ■ a) } 
tr |<t h + tr(Vp) + 8h | - -tr |<r crV(m ■ a) j . 



-tr 



tr(V) 



1 

--tr 
4 



a a tr{Va b 



(A6) 



and m- J m = Y^l b=i ^ abmami >- Varying £ with re- 
spect to m under the condition |m| = 1 one obtains 



m x (B- J -m) = 



(A7) 



This is a "magnetic form" of the equations (|A1[ ). Of 
course i t co uld also be obtained by the direct substitu- 
tion of (A4). The "magnetic field" B contains the infor- 
mation about the random part h + 5h of the HF Hamilto- 

nian whereas the matrix J is the projection on the Pauli 
matrices of the interaction matrix Vijkl taken in the 
subspace of the two states ipi and ■02- For our general 

discussion we can diagonalize J . This will reshuffle the 
components of B which are random anyway. 



The distance Ae between the eigenvalues of (Al) is 
Ae = tr [(020^ - (t)i(j)i)(h + Sh) HF ] where due to or- 
thogonality <p2 $2 = (1 — m - <t)/2. Therefore 



Ae(B, J) = 2 B J m(B, J) 



(A8) 



with m(B, J) - the solution of (hj 

The equations are especially simple in the case when 
the Hamiltonian and the wavefunctions are real. This 
corresponds to the orthogonal ensemble in the terminol- 
ogy of RMT. Only the projections on o\ and (73 remain 
in this case. The HF energy and the HF equation are 
respectively 

£ = B x mi + B 3 ni 3 - ^(Jimj + J 3 m 2 3 ) , (A9) 



TO1-B3 — m 3 Bi + (Ji — J 3 )mim 3 = 



(A10) 



with the constraint 



= 1. In the above J\ and J 3 



are the eigenvalues of J. The solution of (A10) are roots 
of a fourth order equation. Real roots are the extrema of 
the curve defined by the intersection of the surface ( |A9| ) 
and the constraint cylinder m\ + m 3 = 1. Only one real 
root must be selected - that which minimizes £. 
Let us examine the energy spacing Ae 



Ae = 2^J{B l - J^m) 2 + {B 3 - J 3 m 3 y 



(AH) 



and determine when and how level crossings occur. In 
the absence of the interaction (Ji = J 3 = 0) Ae = 
at B\ = B3 = 0, reproducing the standard result. In 
the presence of the interaction the level crossing Ae = 
occurs when 



m 1 =B 1 /J 1 , m 3 = B 3 /J 3 



(A12) 



Due to the normalization constraint this can be fulfilled 
only on an ellipse 



JI + JI 



1 . 



(A13) 



It is therefore sufficient to examine the minimum e nergy 
solution on the ellipse. We note that the conditions (A12) 
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are consistent with equation (A10), i.e. level crossings 
must occur at least for one of the solutions. However 
there is no guarantee that this is also the minimum en- 
ergy solution. In fact a detailed analysis reveals that 
the occurrence of level crossing for the minimum energy 
solution depends on the signs of the eigenvalues J\ and 
J3 as follows: (i) when both signs are positive there are 
four real solutions (two minima and two maxima) on the 
ellipse (A13) provided J3 > 2Ji or Ji > 2J 3 . Other- 
wise there are only two real solutions (one maximum and 
one minimum). The solution with Ae = is always the 
(highest) maximum, (m) when one of the eigenvalues is 
negative and another positive there are four real solu- 
tions on the ellipse (two minima and two maxima). The 
solution with Ae = is either a metastable minimum or 
a maximum. (Hi) only when both eigenvalues are nega- 
tive the solution with Ae = coincides with the absolute 
minimum of the HF energy. 

In the first two cases the minimum energy solution al- 
ways has a non-zero level spacing. Consequently for such 
J\ and J3 the two HF levels never cross and one expects a 
gap in the probability distribution of Ae. In the last case 
of two positive eigenvalues the level spacing distribution 
in the vicinity of Ae = will depen d on the density of 
Ae = const lines close to the ellipse ( A13| ). 

The interaction V which determines Ji and J3 is just 
V A when the two levels e\ and £1 are not degenerate. In 



this case we find from ( ]A6|) J ab = [(VW&- Vi 2 2i)/2]S ab . 
Thus for a positive (negative) definiteEj, i.e., repulsive 
(attractive) interaction Ji = J 3 > (J\ — J3 < 0). 
When the levels are degenerate one must extend the 
two levels treatment to include all degenerate wavefunc- 
tions. However, if the degeneracy is due to spin and 
the wavefunctions are separable spin-space products, one 
can define an equivalent two level problem for the or- 
bital parts with the interaction V — V A + V. For 
this case J ± = (V212 - 2V1221 - Vn22)/2 J«id J 3 = 
[2V1212 - V1221 - (Van + V 22 22)]/2. For such V the signs 
of J\ and J3 are not uniquely related to the nature of the 
interaction. In particular they may both be positive even 
for a repulsive interaction. However, if this happens, we 
expect to find, due to a general argumentEH, another HF 
solution (e.g. lacking spin degeneracy) with lower total 
energy for which the levels do not cross. 

The case J\ — J3 — J is particularly simple since the 
ellipse degenerates into a circle and one obtains 



Ae(B,J) = 2\B + J\ 



(A14) 



where B = \J B\ + i?| . For positive J this result leads, 
after choosing for simplicity P(B) of a Gaussian form, to 
the shifted Wigner distribution, Eq. p2]). For negative 
J (attractive interaction) we find 



J 



P(s;J) = 



2A 2 



(2\J\-s)e 



2A 



2jJl-s \ 2 

+ (2|J| + s)e~ 



This distribution is a sum of two, mirror reflected and 
shifted Wigner distributions resulting from the two roots 
of s = A e(i?, J). The most important feature to note 
in (A15) is the non-zero probability P(s = 0). This is 



a consequence of the linear dependence of Ae(B, J) at 
B=\J\. 
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